{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further Hypothesis Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select this cell and type Ctrl-Enter to execute the code below.\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data = pd.read_csv(\"stars.csv\")\n",
    "type_key = ['Brown Dwarf', 'Red Dwarf', 'White Dwarf', 'Main Sequence', 'Supergiant','Hypergiant']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8 - Goodness of fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your colleague Althea is not a fan of Prof. Xu's temperature scheme. She has her own ideas about how the star classification should be revised.\n",
    "\n",
    "According to Althea's new theory of stellar evolution, there are two subtypes of hypergiant star, with very different temperature distributions. \n",
    "\n",
    "She says that our galaxy has approximately equal numbers of each subtype. \n",
    "\n",
    "If this theory is true, hypergiant temperatures should be distributed according to:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "P(3000K \\leq \\text{T} \\lt 4000K) = 0.5 \\\\\n",
    "P(4000K \\leq \\text{T} \\lt 40000K) = 0.5\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with uniform temperature distributions within each of the two subtypes.\n",
    "\n",
    "Althea asks you to check whether her theory agrees with your data set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question: do hypergiants in the observed data set fit Althea's theory?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pdf associated with the theory is piecewise uniform. We can construct this as a bespoke continuous probability distribution with `scipy.stats`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def althea_F(x):\n",
    "    if(x<3000):\n",
    "        return 0\n",
    "    elif(x<4000):\n",
    "        return 0.5*(x-3000)/(4000-3000)\n",
    "    elif(x<40000):\n",
    "        return 0.5 + 0.5*(x-4000)/(40000-4000)\n",
    "    else:\n",
    "        return 1\n",
    "\n",
    "class althea(stats.rv_continuous):\n",
    "    \"Althea's distribution\"\n",
    "    def _cdf(self, x):\n",
    "        return np.vectorize(althea_F)(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "al = althea()\n",
    "x = np.arange(0,50000,100)\n",
    "\n",
    "y = []\n",
    "for z in x:\n",
    "    y.append(al.pdf(z))\n",
    "\n",
    "xlab = 'temperature'\n",
    "ylab = 'probability density'    \n",
    "\n",
    "ax = plt.axes()\n",
    "ax.set_xlabel(xlab)\n",
    "ax.set_ylabel(ylab)\n",
    "plt.plot(x, y, color='tab:orange', label='theoretical pdf')\n",
    "ax.legend(loc='upper right')\n",
    "plt.show()    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first glance, the histogram for the hypergiants (type 5) does appear to have a similar shape to this pdf:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observed = data[data.type == 5].temperature\n",
    "n = len(observed)\n",
    "\n",
    "xlab = 'temperature'\n",
    "ylab = 'freq'\n",
    "\n",
    "bins = np.linspace(0,50000)\n",
    "ax = plt.axes()\n",
    "ax.set_xlabel(xlab)\n",
    "ax.set_ylabel(ylab)\n",
    "obs_hist = plt.hist(observed, bins, alpha=0.5, label='observed data')\n",
    "ax.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see more clealy, we can overlay a random sample of the same size ($n$=40), drawn from the theoretical distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expected = al.rvs(size=n)\n",
    "\n",
    "xlab = 'temperature'\n",
    "ylab = 'freq'\n",
    "\n",
    "bins = np.linspace(0,50000)\n",
    "ax = plt.axes()\n",
    "ax.set_xlabel(xlab)\n",
    "ax.set_ylabel(ylab)\n",
    "obs_hist = plt.hist(observed, bins, alpha=0.5, label='observed data')\n",
    "exp_hist = plt.hist(expected, bins, alpha=0.5, label='sampled from theoretical distribution')\n",
    "ax.legend(loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How can we test more formally whether the observed data appear to be drawn from the theoretical distribution?\n",
    "\n",
    "Here we have a theoretical distribution *which is not directly related to any of the standard statistical distributions*, so parametric methods such as the Shapiro-Wilk test are not applicable.\n",
    "\n",
    "However, we can still perform a *non-parametric* goodness-of-fit test, by comparing the theoretical cdf with the empirical one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist = np.histogram(observed, bins=100)\n",
    "empirical_dist = stats.rv_histogram(hist)\n",
    "\n",
    "x = np.arange(0,50000,100)\n",
    "\n",
    "y_theoretical = []\n",
    "y_empirical = []\n",
    "\n",
    "for z in x:\n",
    "    y_theoretical.append(al.cdf(z))\n",
    "    y_empirical.append(empirical_dist.cdf(z))\n",
    "\n",
    "xlab = 'temperature'\n",
    "ylab = 'cumulative probability'\n",
    "    \n",
    "ax = plt.axes()\n",
    "ax.set_xlabel(xlab)\n",
    "ax.set_ylabel(ylab)\n",
    "plt.plot(x, y_empirical, label='empirical cdf')\n",
    "plt.plot(x, y_theoretical, label='theoretical cdf')\n",
    "plt.legend(loc='lower right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kolmogorov-Smirnov test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Theory\n",
    "\n",
    "The [*Kolmogorov-Smirnov test*](https://en.wikipedia.org/wiki/Kolmogorov–Smirnov_test) (or K-S test) examines the deviation of the empirical cdf, $F_n(x)$, from the theoretical one, $F(x)$, to assess the goodness-of-fit. The test statistic is the quantity\n",
    "\n",
    "$$\n",
    "D_n= \\sup_x |F_n(x)-F(x)|,\n",
    "$$\n",
    "\n",
    "i.e. the greatest vertical distance between the two curves.\n",
    "\n",
    "The distribution of $D_n$ under the null hypothesis $H_0: F_n(x) = F(x)$ is called the *Kolmogorov distribution*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Application\n",
    "\n",
    "In our example,\n",
    "\n",
    "$H_0$: $F_n(x) = F(x)$  : The observed temperature distribution of hypergiants is described by Althea's theory.\n",
    "\n",
    "$H_1$: $F_n(x) \\ne F(x)$  : The observed temperature distribution of hypergiants is not described by Althea's theory.\n",
    "\n",
    "$\\alpha = 0.05$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Graphically, we can plot $D_n$ for each $x$-value in the plot above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xlab = 'temperature'\n",
    "ylab = 'empirical cdf - theoretical cdf'\n",
    "\n",
    "y_values = np.array(y_empirical)-np.array(y_theoretical)\n",
    "    \n",
    "ax = plt.axes()\n",
    "ax.set_xlabel(xlab)\n",
    "ax.set_ylabel(ylab)\n",
    "plt.plot(x, y_values, color='black')\n",
    "plt.show()\n",
    "\n",
    "D_n = np.max(np.abs(y_values))\n",
    "print('D_n:', D_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test itself is easy in `scipy.stats`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = stats.kstest(observed, al.cdf)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting $p > \\alpha$, so we accept the null hypothesis: the observed data appear to be compatible with Althea's theory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other applications of the K-S test\n",
    "\n",
    "#### Goodness of fit\n",
    "\n",
    "The K-S test can be used as an alternative method for testing normality, or as a goodness-of-fit test for any other theoretical distribution. However, the standard test is only valid if the parameters (e.g. mean and variance) are *not* estimated from the data. \n",
    "\n",
    "If parameters *are* estimated from the data, we will need to use simulation to find an empirical distribution for $D_n$ under $H_0$.\n",
    "\n",
    "\n",
    "#### Two sample test\n",
    "\n",
    "The K-S test is often used to compare two samples (in the absence of a theoretical cdf) to test whether they follow the same distribution. \n",
    "\n",
    "In `scipy.stats`, this is done with the function `ks_2samp`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
